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ON  THE  ADAPTIVE  IMPLEMENTATION  OF  PISARENKO'S  HARMONIC  RETRIEVAL  METHOD* 
V.  U.  Reddy*,  B.  Egardt  and  T.  Kailath* 


XFOSR-TR-  8  2-0476 


SUMMARY 


The  development  of  adaptive  techniques  for  estimating  the  parameters 
of  sinusoidal  signals  in  additive  noise  is  important  in  many  applications. 

The  so-called  adaptive  line  enhancer(ALE)  proposed  by  Widrow  et  al.[l]  has 
been  a  popular  solution.  The  ALE  is  a  tapped-del ay-line  filter  of  some 
fixed  length,  whose  tap  gains  are  recursively  adjusted  by  using  the  so-called 
LMS[1]  algorithm  so  that  they  converge  to  the  solution  of  the  normal  equ¬ 
ations  for  the  one-step  minimum  mean-square-error  prediction  problem.  Another 
popular  solution  uses  a  so-called  ladder  or  lattice  filter  whose  parameters 
are  adjusted  by  using  a  technique,  due  to  Burg[2],  based  on  minimizing  the 
sum  of  certain  forward  and  backward  one-step  prediction  residuals.  Burg’s 
technique  is  an  off-line  one.  Griffiths[3]  merged  the  above  approaches  by 
proposing  a  lattice  filter  whose  coefficients  were  adapted  by  using  the 
LMS  algorithm,  leading  to  what  is  often  called  a  gradient  lattice(or  ladder) 
filter. 

The  above  three  approaches  all  yield  spectral  estimates  with  fairly 
sharp  peaks  but  the  estimates  of  the  sinusoidal  frequencies  invariably 
appear  to  be  biased  when  the  sinusoids  are  observed  in  the  presence  of 
additive  white  noise. 

In  an  attempt  to  improve  the  above  methods  with  respect  to  bias, 

Ulrych  and  Clayton[4]  hav*  proposed  a  least-squares  fitting  of  an  auto¬ 
regressive  model,  based  on  a  criterion  involving  both  forward  and  back¬ 
ward  prediction  errors  but,  unlike  Burgas  method,  without  using  a  ladder 
filter  model.  With  the  help  of  simulation  results  they  have  demonstrated 
that  the  bias  in  the  spectral  estimates  can  be  reduced  significantly 
compared  to  the  Burg  technique.  It  should  be  noted,  though,  that  this  is 
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true  only  for  short  data  lengths,  since  all  the  above  techniques  give 
identical  results  for  large  data  lengths. 

Possibilities  for  improvement  of  the  above  approaches  becomes  apparent 
when  one  considers  that  they  were  all  designed  to  converge  to  the  optimum 
linear  least-squares  solution  for  the  prediction  of  any  random  process, 
and  do  not  specifically  exploit  the  fact  that  the  signal^  are 

sinusoidal.  Pisarenko[5]  was  perhaps  the  first  to  attempt  to  do  this  in  his 
so-called  "harmonic  retrieval"  method,  which  involves  determining  the  mini¬ 
mum  eigenvalue  and  the  corresponding  eigenvector  of  the  covariance  matrix 
of  the  observed  random  process.  Thompson[6]  noted  that  this  eigenvalue- 
eigenvector  computation  was  equivalent  to  a  certain  constrained  gradient- 
search  procedure  for  obtaining  an  adaptive  version  of  Pisarenko's  method. 
Thompson's  simulations  verified  that  the  frequency  estimates  provided  by 
his  procedure  were  unbiased.  However,  the  main  cost  of  this  technique  was 
that  the  initial  convergence  rate  may  be  very  slow  for  certain  poor 
"initial  conditions". 

One  of  the  goals  of  this  paper  is  to  consider  a  way  of  providing  faster 
initial  convergence  by  using  a  different  algorithm  for  the  above  problem. 
Restating  the  constrained  minimization  as  an  unconstrained  nonlinear  pro¬ 
blem,  we  derived  the  following  two  adaptive  algorithms. 

Consider  the  adaptive  filter  with  constrained  coefficients,  as  sugg¬ 
ested  by  Thompson,  shown  in  Fig.l.  The  observed  process,  consisting  of  a 
sum  of  sinusoids  and  white  noise,  is  denoted  by  a  time  series  x(k).  The 
filter  output  e(k)  can  be  expressed  as  the  inner  product 


e(k)=  ATX ( k ) 

(1) 

where 

E^o****** . *  *l-1^ 

(2a) 

X(k)=  [x(k),x(k-l), . ,x(k-L+l)]T 

(2b) 

and 

A=  A/  1!  A|| 

is  the  constrained  unit-norm  vector.  T  denotes  the  matrix  transpose. 


The  adaptation  criterion  for  the  filter  of  Fig.  1  is 


J=  4E[e2(k)]  (3) 

Expressing  e(k)  in  terms  of  unnormalized  weight  vector  A, 
e(k)=  ATX ( k ) /  (f A  f | 


(4) 


(5) 


the  gradient  estimate  at  k-th  sample  instant  is 

^7  J=  1  [  e(k)X(k)  —  e2(k)  A(k-l)  ] 

|lA(k-l)ll  l|A(k-l)|! 

The  time  update  for  the  normalized  i-th  coefficient  is  then  given  by 

%(!<)=  oC(k)[ai(k-l)-yu[e(k)x(k-i)-e2(k)1i(k-l)]}  (6) 

where  0<(k)=  ||A(k-l)||/  ||  A(k)  |} 


and  is  a  positive  scalar  constant.  Equation  (6)  describes  the  constrained 
LMS  algorithm. 

The  stationary  points  of  the  above  algorithm  are  given  by  the  equation 

Eftx(k)-e(k)S]  e(k ) }  =0 

which,  using  Eq.(l),  can  be  simplified  to  give 
E{x(k)XT(k)}  A=  AT  E -^X(k)XT(k)j  A.A 

Clearly,  every  eigenvector  of  the  covariance  matrix  satisfies  this  equation. 
By  a  somewhat  more  involved  argument  it  can  be  shown  that  only  the  eigen¬ 
vector  corresponding  to  the  minimum  eigenvalue  gives  a  stable  stationary 
point. 

b.  Approximate  Deterministic  Least-Squares  Algorithm 

For  an  alternative  algorithm,  we  choose  the  adaptation  criterion  for 
the  filter  of  Fig.  ^as  the  minimization  of 

V=i^e2(s)  (7) 

with  respect  to  the  unit-norm  vector  A.  Because  of  the  constraint  on  the 
weight  vector,  the  minimization  of  (7)  is  a  nonlinear  problem  and  an  exact 
least-squares  solution  does  not  appear  to  exist.  We,  therefore,  derive  an 
approximate  solution  using  a  Gauss-Newton  type  algorithm. 

Once  again,  expressing  e(s)  as 


e(s)=  A'X(s)/  |1  A  1 1 
we  obtain  £- 

X]\J  = _ 1 _ e(s)[X(s)-  Ae(s)] 
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(8) 

(9) 


where 


if*  (s)=  X(s)  -  Ae(s) 


(10) 


From  (8), (9)  and  (10),  the  Gauss-Newton  algorithm  can  be  obtained  as 


f ol 1 ows :  A  A 

S(k)=  £<(k)[A(k-l)-P(k)W)e(k)] 

P(k-l)  V^(k) ^(k)P(k-l) 

P(k)=  P(k-l) - = - 

l+V>(k)P(k-l)  ^(k) 

t(k)=  X(k )-  A(k-l)e(k) 


(11a) 

(lib) 

(11c) 


where oC(k)  is  a  scalar  constant  whose  value  is  chosen  such  that  the  updated 
weight  vector  has  unit  norm.  In  practical  appl ications ,an  exponential  wei¬ 
ghting  with  the  so-called  forgetting  factor;A;is  applied  to  the  data  so 
as  to  track  the  slowly  varying  parameters  of  the  data.  This  weighting  ref¬ 
lects  in  recursion  (lib)  in  two  ways,  i)  The  right-hand-side  is  divided  by 
>v  ,  and  ii)  unity  in  the  denominator  is  replaced  by  X  . 

Simulations  have  been  performed  to  study  the  properties  of  the  above 
two  schemes.  Figures  2  and  3  show  the  spectral  estimates  obtained  with  the 
two  techniques.  Two  sinusoids  of  normalized  frequencies  0.15  and  0.20  are 
used  in  the  examples.  The  most  important  conclusions  that  can  be  drawn  from 
the  results  are  the  following. 

The  least-squares-type  algorithm  has  faster  initial  convergence.  For 
poor  signal -to-noise  ratio,  both  algorithms  perform  similarly  close  to  the 
true  parameters.  The  removal  of  the  bias  in  the  frequency  estimates  is 
slower  than  the  resolution  between  different  frequencies. 
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(b)  Approximate  least-squares  algorithm  (X=0.98,  P(0)=100) 
Fi9-  2  Simulated  spectral  estimates  (SNR=0  dB,  L=7) 
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(b)  Approximate  least-squares  algorithm  (A=0.98,  P(0)=100) 
Fig.  3  Simulated  spectral  estimates  (SNR=12  dB,  L=5) 


